home *** CD-ROM | disk | FTP | other *** search
/ CD Concept 6 / CD Concept 06.iso / mac / UTILITAIRE / RLaB / toolbox / lagrange.r < prev    next >
Text File  |  1994-07-19  |  2KB  |  72 lines
  1. //-------------------------------------------------------------------//
  2. //LAGRANGE  Lagrange interpolation of arbitrary order. Y = LAGRANGE(TAB,X0,N)
  3. //          returns an N-th order interpolated value from table TAB, looking
  4. //          up X0 in the first column of TAB.
  5.  
  6. //          NOTE: TAB's 1st column is checked for monotonicity. It is an
  7. //          error to request a value outside the range of the first column
  8. //          of TAB for X0.
  9.  
  10. //          If the last argument SKIPM is specified, the monotonicity
  11. //          check will be skipped.
  12.  
  13. //          Original Author: Michael F. Saucier 10-16-87
  14. //-------------------------------------------------------------------//
  15.  
  16. lagrange = function (tab, x0, N, skipm)
  17. {
  18.   local (dx, i, j, jj, jmin, lden, lnum, ...
  19.          m, n, seq, sig, tab2, y0);
  20.  
  21.   if (!exist (N)) 
  22.   { 
  23.     error("Wrong number of input arguments."); 
  24.   }
  25.  
  26.   if (!exist (skipm))
  27.   {
  28.     dx = diff (tab[;1]);
  29.     sig = sign (dx[1]);
  30.     if (any (sign (dx) - sig))
  31.     {
  32.       error("First column of the table must be monotonic.");
  33.     }
  34.   }
  35.  
  36.   // Check range of 1st column versus x0
  37.  
  38.   if (tab[1;1] > x0 || tab[tab.nr;1] < x0) 
  39.   {
  40.     error ("lagrange: x0 must be within range of tab");
  41.   }
  42.  
  43.   i = find (tab[;1] == x0);
  44.   if (i != 0) 
  45.   {
  46.     return tab[i;2];
  47.   }
  48.  
  49.   m = tab.nr; n = tab.nc;
  50.   jmin = min(max(min(find(tab[;1] > x0))-fix((N+1)/2),1),m-N);
  51.   tab2 = tab[jmin:jmin+N;];
  52.   jj = 1:N+1;
  53.  
  54.   seq = x0*ones (1, N+1) - tab2[jj;1]';
  55.   lnum = prod (seq) ./ seq;
  56.  
  57.   lden = ones (1, N+1);
  58.   for (i in jj)
  59.   {
  60.     for (j in jj)
  61.     {
  62.       if (j != i) 
  63.       {
  64.         lden[i] = lden[i] * (tab2[i;1] - tab2[j;1]);
  65.       }
  66.     }
  67.   }
  68.   
  69.   y0 = sum (lnum' ./ lden' .* tab2[jj;2]);
  70.   return y0;
  71. };
  72.